      PROGRAM cnorm

C***********************************************************************
C
C     This program computes cumulative normal and its complement
C     interactively.  It handles large and small values of the
C     argument.

C     .. Parameters ..
      INTEGER mxlist
      PARAMETER (mxlist=1001)
      INTEGER strlen
      PARAMETER (strlen=26)
      CHARACTER*(100) msgfmt
      PARAMETER (msgfmt=
     +'('' You are entering a list of arguments for the cumulative norma
     +l'')')
      CHARACTER*(strlen) onest
      PARAMETER (onest='      1.00000000000000')
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ccum,cum,lcum
      INTEGER iarg,iwhich,lstr,narg
      LOGICAL qneg
      CHARACTER ccumst* (strlen),cumst* (strlen)
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION args(mxlist)
C     ..
C     .. External Functions ..
      DOUBLE PRECISION dlanor
      EXTERNAL dlanor
C     ..
C     .. External Subroutines ..
      EXTERNAL cumnor,fmtexp,umdlsd
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs
C     ..


      WRITE (*,FMT=9000)

 9000 FORMAT (' The maximum absolute value of the argument for which ',
     +       /,
     +       ' the cumulative normal can be calculated by this program',
     +       /,' is .678614E8 = 67861400')

   10 WRITE (*,FMT=9010)

 9010 FORMAT (/,/,' Enter choice of action: ',/,5X,
     +       '(0) Exit this program',/,5X,
     +       '(1) Enter single argument for the cumulative normal',/,5X,
     +       '(2) Enter a list of such arguments',/)

          WRITE (*,FMT=*)
     +      'Enter single argument for the cumulative normal'

          narg = 1

          READ (*,FMT=*) args(narg)


      WRITE (*,FMT=9020)

 9020 FORMAT (/,/,T5,'Argument',T29,'Cumulative Normal',T55,
     +       '1 - Cumulative Normal',/,/)


              CALL cumnor(args(narg),cum,ccum)
              WRITE (*,FMT=9030) "Here I am", args(narg),cum,ccum

 9030         FORMAT (A10, 1X,G19.12,1X,G26.15,1X,G26.15)


	STOP
	END



      SUBROUTINE cumnor(arg,result,ccum)
C***********************************************************************
C
C     SUBROUINE CUMNOR(X,RESULT,CCUM)
C
C                               Function
C
C     Computes the cumulative  of    the  normal   distribution,   i.e.,
C     the integral from -infinity to x of
C          (1/sqrt(2*pi)) exp(-u*u/2) du
C
C     X --> Upper limit of integration.
C                                        X is DOUBLE PRECISION
C
C     RESULT <-- Cumulative normal distribution.
C                                        RESULT is DOUBLE PRECISION
C
C     CCUM <-- Compliment of Cumulative normal distribution.
C                                        CCUM is DOUBLE PRECISION
C
C
C     Renaming of function ANORM from:
C
C     Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
C     Package of Special Function Routines and Test Drivers"
C     acm Transactions on Mathematical Software. 19, 22-32.
C
C     with slight modifications to return ccum and to deal with
C     machine constants.
C
C***********************************************************************
C
C
C Original Comments:
C------------------------------------------------------------------
C
C This function evaluates the normal distribution function:
C
C                              / x
C                     1       |       -t*t/2
C          P(x) = ----------- |      e       dt
C                 sqrt(2 pi)  |
C                             /-oo
C
C   The main computation evaluates near-minimax approximations
C   derived from those in "Rational Chebyshev approximations for
C   the error function" by W. J. Cody, Math. Comp., 1969, 631-637.
C   This transportable program uses rational functions that
C   theoretically approximate the normal distribution function to
C   at least 18 significant decimal digits.  The accuracy achieved
C   depends on the arithmetic system, the compiler, the intrinsic
C   functions, and proper selection of the machine-dependent
C   constants.
C
C*******************************************************************
C*******************************************************************
C
C Explanation of machine-dependent constants.
C
C   MIN   = smallest machine representable number.
C
C   EPS   = argument below which anorm(x) may be represented by
C           0.5  and above which  x*x  will not underflow.
C           A conservative value is the largest machine number X
C           such that   1.0 + X = 1.0   to machine precision.
C*******************************************************************
C*******************************************************************
C
C Error returns
C
C  The program returns  ANORM = 0     for  ARG .LE. XLOW.
C
C
C Intrinsic functions required are:
C
C     ABS, AINT, EXP
C
C
C  Author: W. J. Cody
C          Mathematics and Computer Science Division
C          Argonne National Laboratory
C          Argonne, IL 60439
C
C  Latest modification: March 15, 1992
C
C------------------------------------------------------------------
C------------------------------------------------------------------
C  External Function
C------------------------------------------------------------------
C     .. Scalar Arguments ..
      DOUBLE PRECISION arg,ccum,result
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION del,eps,half,min,one,root32,sixten,sqrpi,temp,
     +                 thrsh,x,xden,xnum,xsq,y,zero
      INTEGER i
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION a(5),b(4),c(9),d(8),p(6),q(5)
C     ..
C     .. External Functions ..
      DOUBLE PRECISION d1mach
      EXTERNAL d1mach
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,aint,exp
C     ..
C     .. Data statements ..
C------------------------------------------------------------------
C  Mathematical constants
C
C  SQRPI = 1 / sqrt(2*pi), ROOT32 = sqrt(32), and
C  THRSH is the argument for which anorm = 0.75.
C------------------------------------------------------------------
C------------------------------------------------------------------
C  Coefficients for approximation in first interval
C------------------------------------------------------------------
C------------------------------------------------------------------
C  Coefficients for approximation in second interval
C------------------------------------------------------------------
C------------------------------------------------------------------
C  Coefficients for approximation in third interval
C------------------------------------------------------------------
      DATA one,half,zero,sixten/1.0D0,0.5D0,0.0D0,1.60D0/,
     +     sqrpi/3.9894228040143267794D-1/,thrsh/0.66291D0/,
     +     root32/5.656854248D0/
      DATA a/2.2352520354606839287D00,1.6102823106855587881D02,
     +     1.0676894854603709582D03,1.8154981253343561249D04,
     +     6.5682337918207449113D-2/
      DATA b/4.7202581904688241870D01,9.7609855173777669322D02,
     +     1.0260932208618978205D04,4.5507789335026729956D04/
      DATA c/3.9894151208813466764D-1,8.8831497943883759412D00,
     +     9.3506656132177855979D01,5.9727027639480026226D02,
     +     2.4945375852903726711D03,6.8481904505362823326D03,
     +     1.1602651437647350124D04,9.8427148383839780218D03,
     +     1.0765576773720192317D-8/
      DATA d/2.2266688044328115691D01,2.3538790178262499861D02,
     +     1.5193775994075548050D03,6.4855582982667607550D03,
     +     1.8615571640885098091D04,3.4900952721145977266D04,
     +     3.8912003286093271411D04,1.9685429676859990727D04/
      DATA p/2.1589853405795699D-1,1.274011611602473639D-1,
     +     2.2235277870649807D-2,1.421619193227893466D-3,
     +     2.9112874951168792D-5,2.307344176494017303D-2/
      DATA q/1.28426009614491121D00,4.68238212480865118D-1,
     +     6.59881378689285515D-2,3.78239633202758244D-3,
     +     7.29751555083966205D-5/
C     ..
C------------------------------------------------------------------
C  Machine dependent constants
C------------------------------------------------------------------
      eps = d1mach(3)
      min = d1mach(1)
C------------------------------------------------------------------
      x = arg
      y = abs(x)
      IF (y.LE.thrsh) THEN
C------------------------------------------------------------------
C  Evaluate  anorm  for  |X| <= 0.66291
C------------------------------------------------------------------
          xsq = zero
          IF (y.GT.eps) xsq = x*x
          xnum = a(5)*xsq
          xden = xsq
          DO 10 i = 1,3
              xnum = (xnum+a(i))*xsq
              xden = (xden+b(i))*xsq
   10     CONTINUE
          result = x* (xnum+a(4))/ (xden+b(4))
          temp = result
          result = half + temp
          ccum = half - temp
C------------------------------------------------------------------
C  Evaluate  anorm  for 0.66291 <= |X| <= sqrt(32)
C------------------------------------------------------------------
      ELSE IF (y.LE.root32) THEN
          xnum = c(9)*y
          xden = y
          DO 20 i = 1,7
              xnum = (xnum+c(i))*y
              xden = (xden+d(i))*y
   20     CONTINUE
          result = (xnum+c(8))/ (xden+d(8))
          xsq = aint(y*sixten)/sixten
          del = (y-xsq)* (y+xsq)
          result = exp(-xsq*xsq*half)*exp(-del*half)*result
          ccum = one - result
          IF (x.GT.zero) THEN
              temp = result
              result = ccum
              ccum = temp
          END IF
C------------------------------------------------------------------
C  Evaluate  anorm  for |X| > sqrt(32)
C------------------------------------------------------------------
      ELSE
          result = zero
          xsq = one/ (x*x)
          xnum = p(6)*xsq
          xden = xsq
          DO 30 i = 1,4
              xnum = (xnum+p(i))*xsq
              xden = (xden+q(i))*xsq
   30     CONTINUE
          result = xsq* (xnum+p(5))/ (xden+q(5))
          result = (sqrpi-result)/y
          xsq = aint(x*sixten)/sixten
          del = (x-xsq)* (x+xsq)
          result = exp(-xsq*xsq*half)*exp(-del*half)*result
          ccum = one - result
          IF (x.GT.zero) THEN
              temp = result
              result = ccum
              ccum = temp
          END IF

      END IF

      IF (result.LT.min) result = 0.0D0
      IF (ccum.LT.min) ccum = 0.0D0
C------------------------------------------------------------------
C  Fix up for negative argument, erf, etc.
C------------------------------------------------------------------
C----------Last card of ANORM ----------
      END


      DOUBLE PRECISION FUNCTION d1mach(i)
C***********************************************************************
C
C     DOUBLE PRECISION FUNCTION d1mach(i)
C
C                  Double-precision machine constants
C                         (Truncated version)
C
C
C     Argument:
C
C     I      --> Which constant to return.  See below.
C                INTEGER i
C
C
C     Details:
C
C     D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude.
C
C     D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
C
C     D1MACH( 3) = B**(-T), the smallest relative spacing.
C
C     D1MACH( 4) = B**(1-T), the largest relative spacing.
C
C     D1MACH( 5) = LOG10(B).
C
C     This function has been modified from the original to return values
C     which are specified  by ASCII representations,  which we hope will
C     be  portable to  most systems.   These values  are close to  those
C     which are correct for a 64-bit  DEC Alpha, but truncated slightly.
C     The   truncations were  performed  in  order to  prevent over-  or
C     under-flows, or other problems due to slight system variations.
C
C
C     Note:
C
C     The following values were obtained using  the original D1MACH on a
C     64-bit DEC Alpha (using the DATA statements for little-endian IEEE
C     machines):
C
C     DATA dmach(1)/2.225073858507201E-308/
C     DATA dmach(2)/1.797693134862316E+308/
C     DATA dmach(3)/1.110223024625157E-016/
C     DATA dmach(4)/2.220446049250313E-016/
C     DATA dmach(5)/0.301029995663981/
C
C     The values actually used are truncated versions of these.
C
C***********************************************************************
C     .. Scalar Arguments ..
      INTEGER i
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION dmach(5)
C     ..
C     .. Data statements ..
      DATA dmach(1)/2.22507385850721D-308/
      DATA dmach(2)/1.79769313486231D+308/
      DATA dmach(3)/1.11022302462516D-016/
      DATA dmach(4)/2.22044604925031D-016/
      DATA dmach(5)/0.301029995663981D+0/
C     ..

      IF (i.LT.1 .OR. i.GT.5) THEN
          WRITE (*,FMT='('' D1MACH - I OUT OF BOUNDS'',I10)') i
          STOP

      END IF

      d1mach = dmach(i)

      RETURN

      END



      SUBROUTINE fmtexp(x,string,lstr)
C***********************************************************************
C
C     Writes exp(x) to string even when exp(x) cannot be represented
C     The value is written left adjusted with only as many decimal
C     digits as approximately warranted assuming a double precision
C     number has 16 digits.
C     LSTR is the length of the string.
C
C     Modified to right adjust string in a field 30 char long
C
C***********************************************************************


C     .. Parameters ..
      DOUBLE PRECISION dln10
      PARAMETER (dln10=2.3025850929940456840D0)
      INTEGER strlen
      PARAMETER (strlen=26)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION x
      INTEGER lstr
      CHARACTER string* (strlen)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION digexp,digits,expnt,ldig,y
      INTEGER case,i,ndig,nright
      CHARACTER fmt1*8,expont*22,mantis*22,tmp* (strlen)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,aint,int,log10
C     ..
      string = ' '
      tmp = ' '
      fmt1 = '(F20.XX)'
      mantis = ' '
      expont = ' '

      y = x/dln10

C     EXPNT is the decimal exponent of the number

      expnt = aint(y)
      IF (y.LT.0) expnt = expnt - 1

C      write (*,*) 'EXPNT first calculated: ',EXPNT

C     LDIG is the log10 of the mantissa of the number
C     DIGITS is the mantissa of the number

      ldig = y - expnt
      digits = 10.0D0**ldig

C     NDIG is number of significant digits for mantissa

      digexp = log10(abs(expnt))
      ndig = int(16.0D0-digexp)

C     In case there are no significant digits left return all *

      IF (ndig.LT.1) THEN
          string = '******************************'
          lstr = strlen
          RETURN

      END IF

C     Otherwise let's write out digits to MANTIS with full accuracy

      nright = ndig - 1
      WRITE (fmt1(6:7),FMT='(I2)') nright
      IF (fmt1(6:6).EQ.' ') fmt1(6:6) = '0'
      WRITE (mantis,FMT=fmt1) digits
C      WRITE (*,*) 'Mantissa String: *',MANTIS,'*'

C     Now let's write exponent to EXPONT

      WRITE (expont,FMT='(F18.0)') expnt

C      WRITE (*,*) 'Exponent String: *',EXPONT,'*'

C     Now put it all together in STRING
C     First the mantissa

      lstr = 0
      DO 10 i = 1,20
          IF (mantis(i:i).NE.' ') THEN
              lstr = lstr + 1
              tmp(lstr:lstr) = mantis(i:i)
          END IF

   10 CONTINUE

C     Then 'X10^( '
      tmp(lstr+1:lstr+6) = 'X10^( '
      lstr = lstr + 6

C     Then the exponent

      DO 20 i = 1,20
          IF (expont(i:i).NE.' ' .AND. expont(i:i).NE.'.') THEN
              lstr = lstr + 1
              tmp(lstr:lstr) = expont(i:i)
          END IF

   20 CONTINUE

      tmp(lstr+1:lstr+2) = ' )'
      lstr = lstr + 2

C      WRITE (*,*) ' lstr is: ',lstr

C     Right adjust tmp in string

      string(strlen+1-lstr:) = tmp(1:lstr)

      RETURN

      END



      DOUBLE PRECISION FUNCTION dlanor(x)
C***********************************************************************
C
C      DOUBLE PRECISION FUNCTION DLANOR( X )
C           Double precision Logarith of the Asymptotic Normal
C
C                               Function
C
C      Computes the logarithm of the cumulative normal distribution
C      from abs( x ) to infinity for abs( x ) >= 5.
C
C                              Arguments
C
C      X --> Value at which cumulative normal to be evaluated
C                     DOUBLE PRECISION X
C
C                                Method
C
C      23 term expansion of formula 26.2.12 of Abramowitz and Stegun.
C      The relative error at X = 5 is about 0.5E-5.
C
C                                 Note
C
C      ABS(X) must be >= 5 else there is an error stop.
C
C***********************************************************************

C     .. Parameters ..
      DOUBLE PRECISION dlsqpi
      PARAMETER (dlsqpi=0.91893853320467274177D0)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION x
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION approx,correc,xx,xx2
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION coef(12)
C     ..
C     .. External Functions ..
      DOUBLE PRECISION devlpl,dln1px
      EXTERNAL devlpl,dln1px
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,log
C     ..
C     .. Data statements ..
      DATA coef/-1.0D0,3.0D0,-15.0D0,105.0D0,-945.0D0,10395.0D0,
     +     -135135.0D0,2027025.0D0,-34459425.0D0,654729075.0D0,
     +     -13749310575D0,316234143225.0D0/
C     ..

      xx = abs(x)
      IF (xx.LT.5.0D0) STOP ' Argument too small in DLANOR'

      xx2 = xx*xx
      approx = -dlsqpi - 0.5*xx2 - log(xx)

      correc = devlpl(coef,12,1.0D0/xx2)/xx2
      correc = dln1px(correc)

      dlanor = approx + correc

      RETURN

      END








      SUBROUTINE crlist(xl,xh,nint,qlin,xlist)
C**********************************************************************
C
C     SUBROUTINE crlist(xl,xh,nint,qlin,xlist)
C                      (Double precision version)
C
C     Creates list of  entries between XL  and XH.   There are NINT  + 1
C     such  entries total in XLIST.  IF  QLIN is .TRUE., the entries are
C     evenly spaced; if QLIN is  .FALSE., the entries are evenly  spaced
C     on a log scale.
C
C
C     Arguments
C
C     XL     --> Low value for the list.
C                DOUBLE PRECISION xl
C
C     XH     --> High value for the list.
C                DOUBLE PRECISION xh
C
C     NINT   --> Number of intervals between XL and XH.
C                (Thus the number of values is NINT+1)
C                INTEGER nint
C
C     QLIN   --> If .TRUE., spacing between values is linear.
C                If .FALSE., spacing between values is logarithmic.
C                LOGICAL qlin
C
C     XLIST  <-- List of values between XL and XH.  XLIST should have at
C                    least (NINT+1) elements.
C                DOUBLE PRECISION xlist(*)
C
C
C     Note
C
C     For  pretty  values, NINT should divide  XH  - XL (linear case) or
C     Log(XH/XL) (logarithmic case).
C
C**********************************************************************
C     .. Scalar Arguments ..
      DOUBLE PRECISION xh,xl
      INTEGER nint
      LOGICAL qlin
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION xlist(*)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION xint,xlxl
      INTEGER i
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC dble,exp,log
C     ..
      IF (.NOT. (qlin)) GO TO 20
C
C     Linear List
C
      xint = (xh-xl)/dble(nint)
      DO 10 i = 1,nint + 1
          xlist(i) = xl + dble(i-1)*xint
   10 CONTINUE
      GO TO 40
C
C     Log List
C
   20 xint = log(xh/xl)/dble(nint)
      xlxl = log(xl)
      DO 30 i = 1,nint + 1
          xlist(i) = exp(xlxl+dble(i-1)*xint)
   30 CONTINUE
   40 RETURN

      END



      DOUBLE PRECISION FUNCTION devlpl(a,n,x)
C
C**********************************************************************
C
C     DOUBLE PRECISION FUNCTION DEVLPL(A,N,X)
C              Double precision EVALuate a PoLynomial at X
C
C
C                              Function
C
C
C     returns
C          A(1) + A(2)*X + ... + A(N)*X**(N-1)
C
C
C                              Arguments
C
C
C     A --> Array of coefficients of the polynomial.
C                                        A is DOUBLE PRECISION(N)
C
C     N --> Length of A, also degree of polynomial - 1.
C                                        N is INTEGER
C
C     X --> Point at which the polynomial is to be evaluated.
C                                        X is DOUBLE PRECISION
C
C**********************************************************************
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION x
      INTEGER n
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION a(n)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION term
      INTEGER i
C     ..
      term = a(n)
      DO 10 i = n - 1,1,-1
          term = a(i) + term*x
   10 CONTINUE
      devlpl = term
      RETURN

      END
      DOUBLE PRECISION FUNCTION dln1px(a)
C**********************************************************************
C
C     DOUBLE PRECISION FUNCTION DLN1PX(X)
C               Double precision LN(1+X)
C
C
C                              Function
C
C
C     Returns ln(1+x)
C     Note that the obvious code of
C               LOG(1.0+X)
C     won't work for small X because 1.0+X loses accuracy
C
C
C                              Arguments
C
C
C     X --> Value for which ln(1-x) is desired.
C                                        X is DOUBLE PRECISION
C
C
C                              Method
C
C
C     Renames ALNREL from:
C     DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
C     Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
C     Trans. Math.  Softw. 18 (1993), 360-373.
C
C**********************************************************************
C-----------------------------------------------------------------------
C            EVALUATION OF THE FUNCTION LN(1 + A)
C-----------------------------------------------------------------------
C     .. Scalar Arguments ..
      DOUBLE PRECISION a
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION p1,p2,p3,q1,q2,q3,t,t2,w,x
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,dble,dlog
C     ..
C     .. Data statements ..
      DATA p1/-.129418923021993D+01/,p2/.405303492862024D+00/,
     +     p3/-.178874546012214D-01/
      DATA q1/-.162752256355323D+01/,q2/.747811014037616D+00/,
     +     q3/-.845104217945565D-01/
C     ..
C--------------------------
      IF (abs(a).GT.0.375D0) GO TO 10
      t = a/ (a+2.0D0)
      t2 = t*t
      w = (((p3*t2+p2)*t2+p1)*t2+1.0D0)/ (((q3*t2+q2)*t2+q1)*t2+1.0D0)
      dln1px = 2.0D0*t*w
      RETURN
C
   10 x = 1.D0 + dble(a)
      dln1px = dlog(x)
      RETURN

      END
      SUBROUTINE dlpkdp(dlist,ldlist,loelem,hielem)
C**********************************************************************
C
C     SUBROUTINE dlpkdp(dlist,ldlist,loelem,hielem)
C         DeLete elements from a list  and PacK the list, which contains
C         Double Precision numbers.
C
C
C     Function
C
C     This function  will delete one or  more consecutive numbers from a
C     list of numbers and move  those left up  so that the list is again
C     contiguous.
C
C
C     Arguments
C
C     DLIST  <-> List of   double precision numbers from  which elements
C                    are to be deleted.
C                DOUBLE PRECISION dlist(*)
C
C     LDLIST <-> On input, the original length of DLIST.
C                On output, the   length  of DLIST  with   DLIST(LOELEM)
C                    through DLIST(HIELEM) removed.
C                INTEGER ldlist
C
C     LOELEM --> The index of the first element to be deleted.
C                INTEGER loelem
C
C     HIELEM --> The index of the last element to be deleted.
C                INTEGER hielem
C
C**********************************************************************
C     .. Scalar Arguments ..
      INTEGER hielem,ldlist,loelem
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION dlist(*)
C     ..
C     .. Local Scalars ..
      INTEGER i,idiff
C     ..
      IF (loelem.GT.hielem) STOP ' LOELEM > HIELEM in DLPKDP'
      IF (hielem.GT.ldlist) STOP ' HIELEM > LDLIST in DLPKDP'
      idiff = (hielem-loelem) + 1
      DO 10 i = loelem, (ldlist-idiff)
          dlist(i) = dlist(i+idiff)
   10 CONTINUE
      ldlist = ldlist - idiff
      RETURN

      END





      INTEGER FUNCTION intgt(fmt,low,hi)




C     .. Scalar Arguments ..
      INTEGER hi,low
      CHARACTER fmt* (*)
C     ..
C     .. Local Scalars ..
      INTEGER ifail,ios,nfail
      LOGICAL qfail
C     ..
      DO 10 ifail = 1,3

          IF (fmt.NE.' ') WRITE (*,FMT=fmt)

          READ (*,FMT=*,iostat=ios) intgt

          IF (ios.NE.0) THEN
              WRITE (*,FMT=*)
     +          'The value entered was not a legal integer.'
              GO TO 10

          END IF

          IF (intgt.LT.low .OR. intgt.GT.hi) THEN

              WRITE (*,FMT=9000) intgt,low,hi

 9000         FORMAT (' The value entered: ',I10,' was not in bounds.',
     +               /,' The low value allowed is: ',I10,/,
     +               ' The high value allowed is: ',I10,/,' Try again.')

              GO TO 10

          END IF

          RETURN

   10 CONTINUE

      STOP 'Abort after three attempts to get an integer  in intgt.'

      END
      SUBROUTINE prlist(unit,x,n)
C     .. Scalar Arguments ..
      INTEGER n,unit
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION x(n)
C     ..
C     .. Local Scalars ..
      INTEGER i
C     ..
      DO 10 i = 1,n
          WRITE (unit,FMT=*) i,x(i)
   10 CONTINUE

      RETURN

      END
      LOGICAL FUNCTION qdfdbl(x,y,irow)
C     .. Scalar Arguments ..
      INTEGER irow
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION x(*),y(*)
C     ..
      IF (irow.GT.0) THEN
          qdfdbl = x(irow) .GT. y(irow)

      ELSE
          qdfdbl = x(-irow) .LE. y(-irow)
      END IF

      RETURN

      END
      SUBROUTINE srmdbl(qgr,a,nrowus,nrowdm,ncol,irow)
C**********************************************************************
C
C     SUBROUTINE SRMDBL(QGR,A,NROWUS,NROWDM,NCOL,IROW)
C
C               SoRt Matrix of DouBLe precision values
C
C
C                              Function
C
C
C          Sorts the columns of array A in ascending order acccording
C     to the collating sequence implied by the user supplied
C     function QGR.
C          Note that each element of A must be less than or equal to
C     128 characters in length for the sort to work.  Otherwise the
C     routine will exit with an error indication.
C
C
C                              Arguments
C
C
C     QGR --> The name of a logical function supplied by the user.
C             It has three arguments, X, Y, and IROW, e.g.,
C                    QXGTY = QGR(X,Y,IROW)
C             where X and Y are columns of array A and must be
C             dimensioned accordingly.  IROW is passed through to QGR
C             unchanged from the value sent to this routine.
C                 QGR should return .TRUE. if X is greater than Y in the
C             collating sequence desired and .FALSE. if X is less than
C             or equal to Y.
C                 QGR must be declared to be EXTERNAL in the calling
C             program.
C                                   QGR is LOGICAL
C
C     A --> The array whose columns are sorted.
C                                   A is DOUBLE PRECISION A(NROWDM,NCOL)
C
C     NROWUS --> The number of rows of A that are actually used.  The
C                sort routine works by permuting columns and only rows
C                1..NROWUS of each column are permuted.  The remaining
C                rows are not referenced by this routine.
C                     If a singly dimensioned array is being sorted
C                then NROWUS should be set to 1.
C                                   NROWUS is INTEGER
C
C     NROWDM --> The number of rows declared in A in the program that
C                dimensions it.
C                     If a singly dimensioned array is being sorted
C                then NROWDM should be set to 1.
C                                   NROWDM is INTEGER
C
C     NCOL --> The number of columns of A to be sorted.
C              Columns 1..NCOL are sorted.  It is irrelevent
C              to this routine what the number of columns declared
C              in A is.
C                                   NCOL is INTEGER
C
C     IROW --> An integer argument passed unchanged to QGR.  The intent
C              is that it be the row on which the sort is performed in
C              cases in which that makes sense.  This allows one version
C              of QGR to be used regardless of the row on which the
C              sort is performed.  Note that IROW can be used for
C              any purpose that the user desires.
C                                   IROW is INTEGER
C
C
C                              Note
C
C
C     The algorithm used is based on the article by Robert Sedgewick,
C     "Implementing Quicksort Programs", CACM, Oct 1978, 21:10,847-857.
C
C     If QGR allows ties, the order of the tied columns is not
C     determined by this routine.
C
C**********************************************************************
C***
C
C     DECLARE ARGUMENTS
C
C***
C
C *** VARIABLES
C
C
C  ** PARAMETERS AND CONSTANTS
C
C   * CONSTANTS
C
C     QFALSE  --  MNEMONICALLY-NAMED LOGICAL CONSTANT
C
C
C     QTRUE   --  MNEMONICALLY-NAMED LOGICAL CONSTANT
C
C
C   * PARAMETERS PECULIAR TO THIS APPLICATION
C
C     CHGSRT  --  THRESHOLD ABOVE WHICH QUICKSORT IS FASTER THAN
C                 INSERTION SORT
C
C
C     STCKMX  --  MAXIMUM NUMBER OF ENTRIES IN STACK
C
C
C
C  ** OTHER SIGNIFICANT VARIABLES
C
C     PART    --  INDEX IN ARRAY OF ELEMENT CHOSEN AS PARTITION
C
C
C
C     STCK    --  STACK TO SAVE THE STARTING AND ENDING POSITIONS OF THE
C                 PARTITIONS
C
C
C     STCKCT  --  CURRENT LENGTH OF STACK
C
C
C     SWAP*   --  INDICES IN ARRAY OF ELEMENTS TO BE SWAPPED
C
C
C
C  ** MISCELLANEOUS VARIABLES
C     THE FOLLOWING VARIABLES ARE USED FOR SUCH PURPOSES AS INDICES AND
C     TEMPORARY STORAGE.  EACH ONE IS INTENDED TO BE SIGNIFICANT ONLY
C     WITHIN A WELL-DEFINED BLOCK OF CODE (ALTHOUGH IT MAY APPEAR IN
C     MORE THAN ONE BLOCK), AND ITS MEANING AND USE SHOULD BE APPARENT
C     FROM ITS NAME AND CONTEXT.
C
C     I       --
C
C     J       --
C
C     L       --
C
C     Q       --
C
C     QDONE   --
C
C     R       --
C
C     ICOL    --
C
C     II      --
C
C     JJ      --
C
C     IXROW   --
C
C
C
C
C *** FUNCTIONS AND SUBROUTINES
C
C
C  ** LIBRARY SUBPROGRAMS
C     DETAILED INFORMATION ON EACH OF THE FOLLOWING SUBPROGRAMS CAN BE
C     FOUND IN THE REFERENCE FOR ITS RESPECTIVE LIBRARY.
C
C   * FORTRAN
C     INFORMATION ON EACH OF THE FOLLOWING CAN BE FOUND IN THE FORTRAN
C     LANGUAGE REFERENCE MANUAL.
C
C     MAX     --  MAXIMUM
C
C
C     MIN     --  MINIMUM
C
C
C     LEN     -- LENGTH OF A CHARACTER STRING
C
C
C
C
C
C     *** QUICKSORT
C     .. Parameters ..
      LOGICAL qfalse
      PARAMETER (qfalse=.FALSE.)
      LOGICAL qtrue
      PARAMETER (qtrue=.TRUE.)
      INTEGER chgsrt
      PARAMETER (chgsrt=10)
      INTEGER stckmx
      PARAMETER (stckmx=50)
C     ..
C     .. Scalar Arguments ..
      INTEGER irow,ncol,nrowdm,nrowus
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION a(nrowdm,ncol)
C     ..
C     .. Function Arguments ..
      LOGICAL qgr
      EXTERNAL qgr
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION temp
      INTEGER i,i99997,icol,ixrow,j,l,part,r,stckct,swap1,swap2
      LOGICAL q,qdone
C     ..
C     .. Local Arrays ..
      INTEGER stck(2,stckmx)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC max,min
C     ..
      stckct = 0
      l = 1
      r = ncol
      qdone = (r-l) .LT. chgsrt
   10 IF (qdone) GO TO 320
C
C
C     *** 'MEDIAN-OF-THREE' MODIFICATION.  SEE 'OVERVIEW' ABOVE
C
      swap1 = (l+r)/2
      swap2 = l + 1
C     SWAP-ELEMENTS
      ASSIGN 20 TO i99997
      GO TO 430
C
   20 IF (.NOT. (qgr(a(1,l+1),a(1,r),irow))) GO TO 40
      swap1 = l + 1
      swap2 = r
C     SWAP-ELEMENTS
      ASSIGN 30 TO i99997
      GO TO 430

   30 CONTINUE
   40 IF (.NOT. (qgr(a(1,l),a(1,r),irow))) GO TO 60
      swap1 = l
      swap2 = r
C     SWAP-ELEMENTS
      ASSIGN 50 TO i99997
      GO TO 430

   50 CONTINUE
   60 IF (.NOT. (qgr(a(1,l+1),a(1,l),irow))) GO TO 80
      swap1 = l + 1
      swap2 = l
C     SWAP-ELEMENTS
      ASSIGN 70 TO i99997
      GO TO 430

   70 CONTINUE
   80 i = l + 1
C
C
C     *** PARTITIONING OF SUBFILE
C
      j = r
      part = l
C
   90 GO TO 110

  100 IF (.NOT. (qgr(a(1,part),a(1,i),irow))) GO TO 120
  110 i = i + 1
      GO TO 100

  120 GO TO 140

  130 IF (.NOT. (qgr(a(1,j),a(1,part),irow))) GO TO 150
  140 j = j - 1
      GO TO 130
C
  150 IF (.NOT. (j.LT.i)) GO TO 160
      GO TO 190

      GO TO 180

  160 swap1 = i
      swap2 = j
C     SWAP-ELEMENTS
      ASSIGN 170 TO i99997
      GO TO 430

  170 CONTINUE
  180 GO TO 90
C
  190 swap1 = l
C
      swap2 = j
C     SWAP-ELEMENTS
      ASSIGN 200 TO i99997
      GO TO 430
C
C
C     *** RECURSION STEP
C
  200 IF (.NOT. (max((j-l), (r-i+1)).LE.chgsrt)) GO TO 230
C     *** DEFER FURTHER SORTING OF BOTH SUBFILES TO INSERTION SORT
C
      IF (.NOT. (stckct.EQ.0)) GO TO 210
      qdone = qtrue
      GO TO 220
C     *** POP STACK AND CONTINUE QUICKSORT
  210 l = stck(1,stckct)
      r = stck(2,stckct)
      stckct = stckct - 1
  220 GO TO 310
C
C     *** CONTINUE QUICKSORT ON AT LEAST ONE SUBFILE
C
  230 IF (.NOT. (min((j-l), (r-i+1)).LE.chgsrt)) GO TO 260
C     *** DEFER SMALL TO INSERTION SORT, CONTINUE QUICKSORT ON LARGE
C
      IF (.NOT. ((j-l).GE. (r-i+1))) GO TO 240
C     *** LEFT SUBFILE IS LARGE ONE
C     L = L
      r = j - 1
      GO TO 250
C     *** RIGHT SUBFILE IS LARGE ONE
  240 l = i
C     R = R
  250 GO TO 300
C
C     *** CONTINUE QUICKSORT ON BOTH SUBFILES
C
  260 IF (.NOT. (stckct.GE.stckmx)) GO TO 270
C     *** STACK IS FULL
      STOP ' STACK OVERFLOW IN SRMDBL'

  270 stckct = stckct + 1
C
C     *** PUSH LARGE SUBFILE ONTO STACK, CONTINUE QUICKSORT WITH SMALL
      IF (.NOT. ((j-l).GE. (r-i+1))) GO TO 280
C     *** LEFT SUBFILE IS LARGE ONE
C
      stck(1,stckct) = l
      stck(2,stckct) = j - 1
      l = i
C     R = R
C
      GO TO 290
C     *** RIGHT SUBFILE IS LARGE ONE
C
  280 stck(1,stckct) = i
      stck(2,stckct) = r
C     L = L
      r = j - 1
C
C
  290 CONTINUE
C
  300 CONTINUE
  310 GO TO 10
C
C
  320 CONTINUE
  330 DO 420 i = (ncol-1),1,-1
C
C
C     *** INSERTION SORT
C
C     *** FOR EACH POSITION PRECEDING THE LAST POSITION ...
C
          IF (.NOT. (qgr(a(1,i),a(1,i+1),irow))) GO TO 410
C     *** A(I) NEEDS TO BE MOVED FURTHER OUT
C
          j = i + 1
          GO TO 350

  340     IF (.NOT. (q)) GO TO 380
C
  350     j = j + 1
C
          IF (.NOT. (j.GT.ncol)) GO TO 360
          q = qfalse
          GO TO 370

  360     q = qgr(a(1,i),a(1,j),irow)
C
  370     GO TO 340

  380     DO 400 ixrow = 1,nrowus
              temp = a(ixrow,i)
              DO 390 icol = i,j - 2
                  a(ixrow,icol) = a(ixrow,icol+1)
  390         CONTINUE
              a(ixrow,j-1) = temp
  400     CONTINUE
C
  410     CONTINUE
  420 CONTINUE
C
      RETURN

C     TO SWAP-ELEMENTS
  430 DO 440 ixrow = 1,nrowus
          temp = a(ixrow,swap1)
          a(ixrow,swap1) = a(ixrow,swap2)
          a(ixrow,swap2) = temp
  440 CONTINUE
      GO TO i99997
C
C
      END
      SUBROUTINE umdlsd(fmtmsg,dlist,ldlist,mxlist)
C**********************************************************************
C
C     SUBROUTINE umdlsd(fmtmsg,dlist,ldlist,mxlist)
C      -U-ser -M-o-D-ifies -L-i-S-t of -D-ouble precision numbers
C
C
C                               Function
C
C     Conducts  a dialog  with the user  to create  or  modify a list of
C     double precision values in  DLIST.  An arbitrary number of actions
C     can be taken  to   make the  list.  Each action   is one  of   the
C     following:
C
C          (1) Enter a user specified number of values.
C
C          (2) Enter low value, high value, and number of intervals.
C              The values added will  be equally spaced between  the low
C              and high  values  specified.   Note  that there   will be
C              number of intervals plus  one  values added because  both
C              the low and high values are on the list.
C
C          (3) Like (2), but values are  equally spaced on a logarithmic
C              scale.
C
C          (4) Print the list.
C
C          (5) Delete consecutive elements of the list.
C
C          (6) Delete a single element of the list.
C
C          (7) Sort  the  list of  values in  ascending order and delete
C              duplicates.
C
C     For  user convenience, if the value  of the  first addition to the
C     list has the same value as the last value already on the list, the
C     redundant value is not added.
C
C
C                              Arguments
C
C     FMTMSG --> A Fortran  format containing a  message identifying the
C                  list  being modified.   This message  is printed each
C                  time the list of actions is given.
C                If FMTMSG is blank, then its contents are ignored.
C                CHARACTER FMTMSG* (*)
C
C     DLIST  <-> The list  of double precision   values as specified  or
C                  modified by the user.
C                DOUBLE PRECISION DLIST(MXLIST)
C
C     LDLIST <-> On input, the  length   of DLIST, i.e.,  the  number of
C                  values already specified. If a new  list of values is
C                  to be created, LDLIST should be 0.
C                On output,  the length of the  list of values specified
C                  by the user.
C                INTEGER LDLIST
C
C     MXLIST --> The maximum allowed length  of DLIST.  Attempts to make
C                  list  longer  than  MXLIST will   result  in an error
C                  message and a refusal to take the action.
C                INTEGER MXLIST
C
C**********************************************************************
C     .. Scalar Arguments ..
      INTEGER ldlist,mxlist
      CHARACTER fmtmsg* (*)
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION dlist(mxlist)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION high,low,zx,zy
      INTEGER hidel,i,i99991,i99995,inun,irow,iwhich,lodel,nint,nlo,nup,
     +        outun
      LOGICAL qok
C     ..
C     .. External Functions ..
      INTEGER intgt
      LOGICAL qdfdbl
      EXTERNAL intgt,qdfdbl
C     ..
C     .. External Subroutines ..
      EXTERNAL crlist,dlpkdp,prlist,srmdbl
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,max
C     ..
C     .. Statement Functions ..
      LOGICAL qeqrl
C     ..
C     .. Statement Function definitions ..
      qeqrl(zx,zy) = abs(zx-zy)/max(abs(zx)+abs(zy),1.0D-20) .LT. 1.0D-5
C     ..
      nup = ldlist
      inun = 5
      outun = 6

   10 IF (fmtmsg.NE.' ') WRITE (outun,FMT=fmtmsg)

      WRITE (outun,FMT=9000) ldlist,mxlist

 9000 FORMAT (' The list currently contains ',I4,' values.',/,
     +       ' The maximum number allowed is ',I4,'.')

      WRITE (*,FMT=*)
      WRITE (outun,FMT=9010)

 9010 FORMAT (' Enter (1) ADD values individually specified ',/,
     +       '       (2) ADD values equally spaced between bounds ',/,
     +       '       (3) ADD values as (2) with logarithmic spacing',/,
     +       '       (4) PRINT list',/,
     +       '       (5) DELETE a single element of list',/,
     +       '       (6) DELETE consecutive elements of list',/,
     +       '       (7) SORT list in ascending order and eliminate ',
     +       'duplicate values',/,'       (8) QUIT ')

      iwhich = intgt(' ',1,8)

      IF ((1).NE. (iwhich)) GO TO 60

      WRITE (outun,FMT=9020)

 9020 FORMAT (' Enter the number of values that will be entered')

      nint = intgt(' ',1,mxlist)
      nint = nint - 1
C     CHECK-LEGAL-NUMBERS
      ASSIGN 20 TO i99995
      GO TO 370

   20 IF (.NOT. (qok)) GO TO 50
      WRITE (outun,FMT=9030) nint + 1

 9030 FORMAT (' Enter ',I3,' values')

      READ (inun,FMT=*) (dlist(i),i=nlo,nup)
      IF (.NOT. (ldlist.GT.0)) GO TO 40
C     ADD-TO-LIST
      ASSIGN 30 TO i99991
      GO TO 340

   30 CONTINUE
   40 ldlist = nup
   50 GO TO 330

   60 IF ((2).NE. (iwhich)) GO TO 100

      WRITE (outun,FMT=9040)

 9040 FORMAT (' Enter low value, high value, number of intervals')

      READ (inun,FMT=*) low,high,nint
C     CHECK-LEGAL-NUMBERS
      ASSIGN 70 TO i99995
      GO TO 370

   70 IF (qok) CALL crlist(low,high,nint,.TRUE.,dlist(nlo))
      IF (.NOT. (ldlist.GT.0)) GO TO 90
C     ADD-TO-LIST
      ASSIGN 80 TO i99991
      GO TO 340

   80 CONTINUE
   90 ldlist = nup
      GO TO 330

  100 IF ((3).NE. (iwhich)) GO TO 140

      WRITE (outun,FMT=9040)
      READ (inun,FMT=*) low,high,nint
C     CHECK-LEGAL-NUMBERS
      ASSIGN 110 TO i99995
      GO TO 370

  110 IF (qok) CALL crlist(low,high,nint,.FALSE.,dlist(nlo))
      IF (.NOT. (ldlist.GT.0)) GO TO 130
C     ADD-TO-LIST
      ASSIGN 120 TO i99991
      GO TO 340

  120 CONTINUE
  130 ldlist = nup
      GO TO 330

  140 IF ((4).NE. (iwhich)) GO TO 160

      IF (.NOT. (ldlist.GT.0)) GO TO 150
      CALL prlist(outun,dlist,ldlist)
  150 GO TO 330

  160 IF ((5).NE. (iwhich)) GO TO 190

      IF (.NOT. (ldlist.GT.0)) GO TO 170
      WRITE (outun,FMT=9050)

 9050 FORMAT (' Enter position  in array for value to be',/,' deleted.')

      lodel = intgt(' ',1,mxlist)
      hidel = lodel
      CALL dlpkdp(dlist,ldlist,lodel,hidel)
      GO TO 180

  170 WRITE (*,FMT=*) 'You must create a list before changing it.'
  180 GO TO 330

  190 IF ((6).NE. (iwhich)) GO TO 220

      IF (.NOT. (ldlist.GT.0)) GO TO 200
      WRITE (outun,FMT=9060)

 9060 FORMAT (' Enter low then high positions in array for',/,
     +       ' values to be deleted.')

      lodel = intgt(' ',1,mxlist)
      hidel = intgt(' ',1,mxlist)
      CALL dlpkdp(dlist,ldlist,lodel,hidel)
      GO TO 210

  200 WRITE (*,FMT=*) 'You must create a list before changing it.'
  210 GO TO 330

  220 IF ((7).NE. (iwhich)) GO TO 290

      IF (.NOT. (ldlist.GT.0)) GO TO 270
      irow = 1
      CALL srmdbl(qdfdbl,dlist,1,1,ldlist,irow)
      i = 1
      GO TO 240

  230 i = i + 1
  240 IF (i.GT.ldlist-1) GO TO 260
      IF (.NOT. (dlist(i).EQ.dlist(i+1))) GO TO 250
      CALL dlpkdp(dlist,ldlist,i+1,i+1)
  250 GO TO 230

  260 GO TO 280

  270 WRITE (*,FMT=*) 'You must create a list before sorting it.'
  280 GO TO 330

  290 IF ((8).NE. (iwhich)) GO TO 320

      IF (.NOT. (ldlist.GT.0)) GO TO 300
      RETURN

      GO TO 310

  300 WRITE (outun,FMT=9070)

 9070 FORMAT (' No values have been placed in the list',/,
     +       ' You can''t quit yet.')

      qok = .FALSE.
  310 GO TO 330

  320 CONTINUE
  330 GO TO 10

C     TO ADD-TO-LIST
C
C
C     Check for equal values at action break point (where old list joins
C     new entries.  If  the values are  equal, eliminate  the second one
C     and adjust NUP accordingly.
C
C
  340 IF (.NOT. (qeqrl(dlist(ldlist),dlist(ldlist+1)))) GO TO 360
      DO 350 i = nlo,nup
          dlist(i-1) = dlist(i)
  350 CONTINUE
      nup = nup - 1
      ldlist = nup
  360 GO TO i99991
C     TO CHECK-LEGAL-NUMBERS

C
C     Check   for number of   intervals.  It must be  at  least one.  If
C     number of intervals is correct, sets the  correct array bounds for
C     the new entries.  Also checks  for exceeding the maximum number of
C     entries.
C

  370 IF (.NOT. (nint.LT.0)) GO TO 380
      WRITE (outun,FMT=9080)

 9080 FORMAT (' Number of intervals must be at least 1')

      GO TO i99995

  380 nlo = ldlist + 1
      nup = ldlist + nint + 1
      qok = nup .LE. mxlist
      IF (.NOT. (qok)) GO TO 390
      GO TO i99995

      GO TO 400

  390 WRITE (outun,FMT=9090)

 9090 FORMAT (' Attempt to add too many values to list',/,' Try again.')

      nup = ldlist
  400 GO TO i99995

      END
